/*
 * ElectrolyteSolution.cpp
 *
 *  Created on: 19 Aug 2011
 *      Author: Allan
 */

#include "ElectrolyteSolution.h"

// GeoReact includes
#include <Geochemistry/Activity/AqueousActivity.h>
#include <Geochemistry/GeoDatabase.h>
#include <Utils/Algorithms.h>
#include <Utils/SpeciesMap.h>

ElectrolyteSolution::ElectrolyteSolution()
{}

ElectrolyteSolution::ElectrolyteSolution(const vector<string>& species) :
species(species)
{
	// Determine the aqueous ions from the given in (species)
	ions = Filter(species, GeoDatabase::IsIonic);
	
	// Determine the aqueous complexes from the given in (species)
	complexes = Filter(species, GeoDatabase::IsComplex);
	
	// The species map object to determine the species indexes
	const SpeciesMap map(species);
	
	// Set the indexes of the ions
	indexes_ions = map[ions];
	
	// Set the indexes of the complexes
	indexes_complexes = map[complexes];
	
	// Set the index of the water species
	iH2O = map["H2O"];
	
	// Initialize the electrical charges of the ions
	zi.resize(ions.size());
	
	for(unsigned i = 0; i < zi.size(); ++i) 
		zi[i] = GeoDatabase::ElectricalCharge(ions[i]);
	
	// Initialize the stoichiometric matrix between the aqueous complexes and ions
	InitializeStoichiometricMatrixComplexesIons();
}

const vector<string>& ElectrolyteSolution::GetSpecies() const
{
	return species;
}

const vector<string>& ElectrolyteSolution::GetIons() const
{
	return ions;
}

const vector<string>& ElectrolyteSolution::GetComplexes() const
{
	return complexes;
}

const VectorXd& ElectrolyteSolution::GetIonicCharges() const
{
	return zi;
}

const Index	ElectrolyteSolution::GetWaterIndex() const
{
	return iH2O;
}

const double ElectrolyteSolution::EffectiveIonicStrength(const VectorXd& n) const
{
	double Ie = 0.5 * (zi.array() * zi.array() * IonicMolalities(n).array()).sum();
	
	return Ie;
}

const double ElectrolyteSolution::StoichiometricIonicStrength(const VectorXd& n) const
{
	double Is = 0.5 * (zi.array() * zi.array() * StoichiometricIonicMolalities(n).array()).sum();
	
	return Is;
}

const VectorXd ElectrolyteSolution::IonicMolalities(const VectorXd& n) const
{
	// The ionic molalities
	VectorXd mi(ions.size());
	
	// The molality factor to trasnform number of moles into molality
	const double molality_factor = 55.508 / n[iH2O]; 
	
	// Loop over all the ions
	for(unsigned i = 0; i < ions.size(); ++i) mi[i] =  molality_factor * n[indexes_ions[i]];
	
	return mi;
}

const VectorXd ElectrolyteSolution::StoichiometricIonicMolalities(const VectorXd& n) const
{
	// Check if there is at least one aqueous complex, otherwise return the ionic molalities
	if(complexes.empty()) return IonicMolalities(n);
	
	// The number of ions (Ni) and complexes (Nc) in the electrolyte solution
	const unsigned Ni = ions.size();
	const unsigned Nc = complexes.size();

	// The moles of the ions (ni)
	VectorXd ni(Ni); for(unsigned i = 0; i < Ni; ++i) ni[i] = n[indexes_ions[i]];
	
	// The moles of the complexes (nc)
	VectorXd nc(Nc); for(unsigned c = 0; c < Nc; ++c) nc[c] = n[indexes_complexes[c]];

	// The molality factor to trasnform number of moles into molality
	const double molality_factor = 55.508 / n[iH2O]; 
	
	// Return the molality factor times the stoichiometric number of moles of each ion
	return molality_factor * (ni + vci.transpose() * nc);
}

const VectorXd ElectrolyteSolution::Molalities(const VectorXd& n) const
{
	// The molality factor to trasnform number of moles into molality
	const double molality_factor = 55.508 / n[iH2O]; 
	
	return molality_factor * n;
}

void ElectrolyteSolution::InitializeStoichiometricMatrixComplexesIons()
{
	// Check if the data members (complexes) and (ions) have been initialized
	assert(!complexes.empty()); assert(!ions.empty());
	
	// Initialise the stoichiometric matrix of the complexes-ions
	vci = MatrixXd::Zero(complexes.size(), ions.size());
	
	// Loop over all the complexes in the electrolyte solution
	for(unsigned c = 0; c < complexes.size(); ++c)
	{
		auto opt_complex_data = GeoDatabase::GetInstance().FindComplex(complexes[c]); assert(opt_complex_data);
		
		// Loop over all the ions that compose the current aqueous complex
		BOOST_FOREACH(auto pair, opt_complex_data->composition)
		{
			const double stoichiometry = pair.first;
			
			const string ion = pair.second;
			
			const unsigned i = IndexOf(ion, ions); // The local index of the current (ion) in the list of (ions)
			
			vci(c, i) = stoichiometry;
		}
	}
}
